package com.lifeonwalden.copula;

import com.lifeonwalden.copula.bean.OptimizeBounds;
import com.lifeonwalden.copula.bean.OptimizeResult;
import com.lifeonwalden.copula.constant.CopulaConstants;

import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.function.Function;

public interface Optimize {
    /**
     * @param func
     * @param bounds
     * @param xatol             Absolute error in solution `xopt` acceptable for convergence.
     * @param maxIterationTimes Maximum number of iterations to perform.
     */
    static OptimizeResult minimizeScalarBounded(Function<BigDecimal, BigDecimal> func, OptimizeBounds bounds, BigDecimal xatol, int maxIterationTimes) {
        BigDecimal lowerBound = bounds.lower, upperBound = bounds.upper;
        if (lowerBound.compareTo(upperBound) > 0) {
            throw new RuntimeException("The lower bound exceeds the upper bound.");
        }

        BigDecimal e = BigDecimal.ZERO;
        BigDecimal rat = BigDecimal.ZERO;

        BigDecimal sqrtEps = BigDecimal.valueOf(Math.sqrt(2.2e-16));
        BigDecimal goldenMean = CopulaConstants.POINT_FIVE.multiply(CopulaConstants.THREE.subtract(BigDecimal.valueOf(Math.sqrt(5.0))));

        BigDecimal fulc = lowerBound.add(goldenMean.multiply(upperBound.subtract(lowerBound)));
        BigDecimal x = fulc, xf = fulc, nfc = fulc;

        BigDecimal fx = func.apply(x);
        BigDecimal ffulc = fx, fnfc = fx;
        BigDecimal fu;

        BigDecimal xm = CopulaConstants.POINT_FIVE.multiply(lowerBound.add(upperBound));

        BigDecimal tol1 = sqrtEps.multiply(xf.abs()).add(xatol.divide(CopulaConstants.THREE, 16, RoundingMode.DOWN));
        BigDecimal tol2 = CopulaConstants.TWO.multiply(tol1);

        int iterationTimes = 1;
        boolean fundSolution = true;
        while (xf.subtract(xm).abs().compareTo(tol2.subtract(CopulaConstants.POINT_FIVE.multiply(upperBound.subtract(lowerBound)))) > 0) {
            BigDecimal golden = BigDecimal.ONE;

            // Check for parabolic fit
            if (e.abs().compareTo(tol1) > 0) {
                golden = BigDecimal.ZERO;
                BigDecimal r = xf.subtract(nfc).multiply(fx.subtract(ffulc));
                BigDecimal q = xf.subtract(fulc).multiply(fx.subtract(fnfc));
                BigDecimal p = xf.subtract(fulc).multiply(q).subtract(xf.subtract(nfc).multiply(r));
                q = CopulaConstants.TWO.multiply(q.subtract(r));

                if (q.compareTo(BigDecimal.ZERO) > 0) {
                    p = p.multiply(CopulaConstants.MINUS_ONE);
                }
                q = q.abs();
                r = e;
                e = rat;

                // Check for acceptability of parabola
                if (p.abs().compareTo(CopulaConstants.POINT_FIVE.multiply(q).multiply(r).abs()) < 0
                        && p.compareTo(q.multiply(lowerBound.subtract(xf))) > 0
                        && p.compareTo(q.multiply(upperBound.subtract(xf))) < 0) {
                    rat = p.divide(q, 16, RoundingMode.DOWN);
                    x = xf.add(rat);

                    if (x.subtract(lowerBound).compareTo(tol2) < 0 || upperBound.subtract(x).compareTo(tol2) < 0) {
                        rat = tol1.multiply(BigDecimal.valueOf(xm.subtract(xf).signum() + (((xm.subtract(xf)).compareTo(BigDecimal.ZERO) == 0) ? 1 : 0)));
                    }
                } else {
                    // do a golden-section step
                    golden = BigDecimal.ONE;
                }
            }

            // do a golden-section step
            if (golden.compareTo(BigDecimal.ZERO) > 0) {
                if (xf.compareTo(xm) >= 0) {
                    e = lowerBound.subtract(xf);
                } else {
                    e = upperBound.subtract(xf);
                }

                rat = goldenMean.multiply(e);
            }

            BigDecimal si = BigDecimal.valueOf(rat.signum() + ((rat.compareTo(BigDecimal.ZERO) == 0) ? 1 : 0));

            x = xf.add(si.multiply(rat.abs().max(tol1)));
            fu = func.apply(x);

            iterationTimes++;

            if (fu.compareTo(fx) <= 0) {
                if (x.compareTo(xf) >= 0) {
                    lowerBound = xf;
                } else {
                    upperBound = xf;
                }

                fulc = nfc;
                ffulc = fnfc;
                nfc = xf;
                fnfc = fx;
                xf = x;
                fx = fu;
            } else {
                if (x.compareTo(xf) < 0) {
                    lowerBound = x;
                } else {
                    upperBound = x;
                }

                if (fu.compareTo(fnfc) <= 0 || nfc.compareTo(xf) == 0) {
                    fulc = nfc;
                    ffulc = fnfc;
                } else if (fu.compareTo(ffulc) <= 0 || fulc.compareTo(xf) == 0 || fulc.compareTo(nfc) == 0) {
                    fulc = x;
                    ffulc = fu;
                }
            }

            xm = CopulaConstants.POINT_FIVE.multiply(lowerBound.add(upperBound));
            tol1 = sqrtEps.multiply(xf.abs()).add(xatol.divide(CopulaConstants.THREE, 16, RoundingMode.DOWN));
            tol2 = CopulaConstants.TWO.multiply(tol1);

            if (iterationTimes > maxIterationTimes) {
                fundSolution = false;
                break;
            }
        }

        return new OptimizeResult().setFun(fx).setX(xf).setNfev(iterationTimes).setNit(iterationTimes).setFundSolution(fundSolution);
    }
}
